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Abstract 

A cluster density matrix is introduced in the form suitable for the self- 
consistent calculation of relevant thermodynamic averages for the Ising model 
with spin S = 1/2. On this basis, derivation of the Gibbs free-energy for the 
effective field theory of Honmura and Kaneyoshi is presented. 
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1 Introduction 



The invention of the differential operator method by Honmura and Kaneyoshi 
(H-K) pQ in 1979, as an approach to the Effective Field Theory (EFT), made 
possible the theoretical investigations of various magnetic models on a large 
scale [2, El HI E], E], [7] . First of all, the method found application in combina- 
tion with the exact Callen-Suzuki statistical identities [HE], giving possibility 
for taking into consideration the correlation function. In its simplest version, 
with only spin- autocorrelations taken into account, the method is equivalent 
to the Zernike [10J, or to 1st Matsudaira [11] approximation. The perspec- 
tives and possible development of the EFT method has been outlined in 
Ref.[12j. In Ref.[T3] the equivalence of the differential operator method to 
another, so-called the integral operator approach, has been discussed. 

Among wide applications of the H-K differential operator method, sev- 
eral important problems should be mentioned. It is, among other things, the 
theory of disordered diluted alloys [31 EH EES] , random field problems [H)J [H] , 
amorphization [TH [T7] IT8] and surface magnetism [HQS]. In the equivalent, 
integral operator version, the method has been applied to studies of thin 
magnetic films [201 IZI] ■ At the same time, great efforts were put into appli- 
cation of H-K method to higher spins problem [221 [23]- In this area, both 
the pure spin models (containing the spins of the same kind), for instance 
the isotropic Blume- Emery-Griffiths model [23], as well as the mixed-spin 
systems [6j [25] have been considered. A lot of attention has also been paid 
to the transversal models [7J [2S1 [27] . 

Despite wide applications and almost 30 year long history of the method, 
the Gibbs free-energy in H-K approximation has not been found so far. Some 
attempts in this direction [28, 29J cannot be considered as satisfactory be- 
cause they fail to lead to the self-consistent thermodynamics. The only ther- 
modynamic potential which has been correctly defined up to now is the in- 
ternal energy [3D] . 

On the other hand, Gibbs free-energy is of particular importance for the 
thermodynamic description. Not least because in case of the 1st order phase 
transitions, or when the several mathematical solutions can occur, one has to 
have a clear criterion to choose up the physical solution corresponding to the 
minimum of free energy. Such problems often occur in connection with the 
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structural disorder in the system, especially in case of higher spins (S > 1/2) 
and the presence of anisotropy. 

The aim of the present paper is a construction of the Gibbs energy in 
H-K approximation, starting from a new form proposed for the cluster den- 
sity matrix. In order to describe the method in detail, the simplest case, 
i.e., the crystalline Ising model with spin S = 1/2, is considered. Further 
progress of the proposed approach (for instance, implementation for higher 
spins) should be possible on this basis. In the derivation of the Gibbs free- 
energy we demand that some known EFT equations for the magnetization 
and NN correlation function remain fulfilled. 

Among the results, the formulas for the entropy and Gibbs free-energy 
in H-K approximation are presented. Some testing numerical calculations 
illustrate the Gibbs free-energy behaviour upon the temperature and, by 
comparison with other mean-field results, serve as a control for the presented 
formalism. 



2 The cluster density matrix 

Below, we consider the crystalline Ising model in an arbitrary dimension 
with spin S = 1/2. We take into account the standard nearest-neighbour 
(NN) exchange interaction and the external field terms. Hamiltonian can be 
written in the form: 

n = -jj:s l s 3 -hj:s 3 (i) 

where Si, Sj = ±1/2 are the ^-components of the spins in i-th and j-th site, 
respectively, J is the exchange interaction for NN, and h = —g^H stands 
for an external magnetic field H oriented in ^-direction. 

In the system we distinguish a (z + l)-atomic cluster, consisting of a 
central spin Si and ^-spins Sj being NN of the central spin. The cluster is 
embedded in an external magnetization field. We restrict our considerations 
for those lattices where the spins Sj in the neighbourhood of Si are not the 
mutual NN. This excludes, for instance, the triangular and FCC lattices. 
However, it should be mentioned that in the case of triangular lattice, the 
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stair-triangle transformation [3TJ E2] can be applied for mapping onto hon- 
eycomb lattice (and vice- versa). 

The density matrix for such a cluster can be proposed in the following 
factorial form: 



Pi{NN} = \ ^ + Si tanh 



jei 



\{[\ + 2mS 3 ) (2) 



jei 



def 

where m = (Si) = (Sj) is a self-consistent parameter being the thermo- 
dynamic mean value of a given spin from the cluster and corresponding to 
the local magnetization. The form of Pi{NN} can be interpreted as an outer 
product of 2 x 2 sigle-site matrices, giving 2 Z+1 x 2 Z+1 diagonal matrix for 
the cluster. It is easy to check that the cluster density matrix (j2J) satisfies 
the normalization condition: 

^i{NN}Pi{NN} = 1 (3) 

where the trace Tr^AW} is performed over the central z-th spin and all the 
NN spins. 

The cluster density matrix (J2J) can be used for calculation of the local 
magnetization m: 

m = Tr i{NN} (SiPi{NN}) (4) 
and the NN correlation function c, defined as: c = (SiSj), namely: 

C = Tli{NN} (SiSjPi{NN}) (5) 

For convenience of calculations of m and c let us note that the matrix pi{NN} 
can be partially reduced. For instance, performing the trace over all NN-spins 
we obtain the single-site density matrix pi 

pi = Tr {NN} p i{NN} = - + 2 (a + a x m + . . . + a z m z ) Si (6) 

where a , ...,a z are the temperature and field dependent coefficients, whose 
explicit form for a given z-number are presented in the Appendix. 
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On the other hand, the single-site density matrix can be generally pre- 
sented in the form (see e.g.Ref.[33j): 

Pl = i + 2mS t (7) 

Hence, by comparison of eqs. (jSJ) and (JTj) we obtain: 

m = ao + aim + . . . + a z m z (8) 

which is the EFT equation for the magnetization in H-K approximation pQ. 
This equation has previously been derived in a different way, from the exact 
Callen-Suzuki identity (see, for instance, in Refs.[T|[T3]) through the operator 
techniques. 

Another partial reduction of Pi{NN} given by eq. ([2]) can be obtained by 
performing the trace over (z — l)-spins from NN, with except of one j-th 
spin. This procedure leads to the pair density matrix p^, namely: 

Pij = T?{NN^j}Pi{NN} = t + (ao + aim + . . . + a z m z ) Si + mSj 

+4 (b + bim + . . . + b z m z ) SiSj (9) 

In (jHJ) bo, b z are the coefficients depending on the ^-number, and are given 
in the Appendix. On the other hand, the pair density matrix should have 
the following symmetrical form [33] : 

Pij — — + mSi + mSj + AcSiSj (10) 

Hence, by comparison of eq. (jHJ) and ( |T0l) we obtain: 

c = b + b x m + . . . + b z m z (11) 

together with eq. (JS]). Eq. ( TTTT) reproduces the EFT result in H-K approxima- 
tion for the NN correlation function, obtained from the exact Callen-Suzuki 
identity (see, for instance, Refs.[13l [30]). Thus, we see that some basic EFT 
equations for m and c which are present in literature can now be alternatively 
re-derived from reducibility of the cluster density matrix Pi{ nn} @. 
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3 The cluster entropy 



Calculation of the entropy presents usually the most difficult task in statisti- 
cal thermodynamics and has not been done yet in the H-K approximation. In 
order to calculate the cluster entropy we will make use of the cluster density 
matrix Pi{nn}- However, in this case the matrix should be transformed first 
to a more convenient, exponential form. For this purpose, we will use the 
exact identity for spin 5 = 1/2: 



' (3A 



pA" 



cosh — — + 2 Si sinh — — 



(12) 



With the help of (Tl2~i) we obtain the next useful relationship: 



— h 2mSi 
2 3 



where 



m 



2 cosh (f) 
(3A\ 



- tanh , 
2 V 2 



(13) 



(14) 



Eq. (Till) defines the A-parameter, which can be unambigously determined 
when the magnetization m is calculated from eq. ((H]). With the help of (fl3l) 
the cluster density matrix (T5]) can be presented in the exponential form: 



pi{NN} 



exp 


P 




St 


2 cosh 


f 







exp PA Sj 



2* cosh 2 [f ) 



(15) 



As a next step, the cluster entropy ai{NN} can be calculated from the formula: 

&i{NN} = (In Pi{NN} ) = — k B Tri{NN} (pi{NN}\n Pi{NN}) (16) 



where Pi{NN} is given by f|T5|) . As a result we obtain: 



&i{NN} = ~j< ( Z J C + hm + zAm) + zk B In 



2 cosh ( ^~ 



+ks ( In < 2 cosh 



(17) 
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The average value of the logarithm function remaining in the last term of 
eq. ffTT|) can be conveniently calculated with the help of Pi{NN} in the form 
of (J2]). The final result is a polynomial of m: 



In < 2 cosh 



jj2s 3 + h 



c + C\m + . . . + c z m z 



(18) 



where c$,...,c z coefficients are given in the Appendix. Finally, the cluster 
entropy ffTTl) is given by the equation: 

' f3A^ 



T(Ti{NN} — —zJc — hm — zAm + zknTln 
+k B T (c + Cim + . . . + c z m z ) 



2 cosh 



(19) 



It is seen from eq. ( I19p that the cluster entropy can be obtained for a given 
structure (^-number) when m and c are calculated first as the solutions of 
EFT equations ([8]) and ffTTl) . respectively. 



4 The Gibbs free-energy 

As stated before, a correct Gibbs potential has not been obtained yet in the 
H-K method. Here we postulate that the Gibbs free-energy of a small cluster 
interacting with its neighbourhood and with an external field can be written 
in the form: 

Gi{nn} = Go — zJc — zmX — (z + 1) hm — Ta^y (20) 

where Go is some constant part of the internal energy ( which doesn't depend 
on the temperature), —zJc is the internal energy of the cluster connected 
with the internal NN corelations only, and A-parameter corresponds to the 
effective-field interaction of the cluster boundary with the rest of the system. 
In eq. (120 p we have also the Zeeman term —{z + l)hm, as well as the entropic 
term — Ta^N}, given by eq. (|T9|) . Then, substituting eq. (|T9|) into (|20|) the 
cluster Gibbs free-energy can be written as: 

Gi{NN} = G — zm(X + A + h) — zk-^T dn 

-k B T (c + dm + . . . + c z m z ) (21) 



2 cosh 
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where A is given by eq. (1141) . The new parameter A can now be determined 
from the condition that the Gibbs energy in equilibrium should be minimized 
with respect to it, or equivalently, G^nn} should be in a minimum with 
respect to the magnetization, since dG^NN^/dX = {dGi^^/drrvj (dm/dX). 
Thus we demand that: 

0. (22) 



9Gi{NN} 



dm 

Eq. (1221) presents the necessary extremum condition for the free-energy, ex- 
ploited in various effective field approaches [281 E3] ■ in our case, from this 
condition we obtain the equation for A, namely: 



dX 
am 



+ X = A - h - -k B T (ci + 2c 2 m + ... + zc z m z -^ (23) 



The above differential equation can be numerically solved only when the 
magnetization m is calculated first from eq. ([8]). The boundary condition 
for A can easily be established for the paramagnetic state. Namely, if we 
put h — 0, m = and A = 0, then from (I2"3"j) we obtain A = 0. Thus, we 
can assume as a boundary condition that A = above the Curie temperature. 



After calculations of A vs. temperature the constant G in eq. (120]) can 
be found from the ground state energy. Namely, we should demand that the 
internal energy per site in the ground state (at T = 0) when calculated from 
( F2H|) is exact, as it results from states of the hamiltonian ([I]). Thus we have 
the condition: 

( G i{N N} \ _ G -(\J+lX )z _ n H )\ = _Jz 

\ z + 1 J t=o z+1 ~ \ N J T =o 8 1 ' 

h=0 h=0 

where Ao = A(T = 0). Solving eq. (124"]) we determine the constant: 

(25) 



Go - 2 



Ao -j{z 



It is worth noticing that the Gibbs energy (121!) is fully described in terms of 
h and T as the primary variables, whereas the intermediate variable m(h, T) 
should be treated solution of eq. (|Sj). 
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It is necessary to check whether the Gibbs energy derivatives properly 
reproduce basic thermodynamic quantities. Differentiating eq. (1211) over h, 
with the help of fl23l) . we obtain after calculation: 



{NN} \ 



dh J, 



-{z + l)m (26) 



which gives the cluster magnetization with m in the form ([8]). During the 
differentiation we have used the following relationship valid for the coefficients 
given in the Appendix: 

(t) T ^» < 27 > 

On the other hand, by differentiation of eq. (1211) over T (with the help of 
) we are able to reproduce the cluster entropy: 



'dG, 



i{NN} 



dT 



-&i{NN} (28) 



where &i{NN} is given by eq. (|T9|) . In this differentiation we made use of 
another useful relationship for the coefficients: 



' dc, 



n 



dT / 1 k-aT 2 



(a n h + zJb n ) (29) 



In a consequence of the first derivatives check, the second derivatives of 
the Gibbs potential (!2~T1) will lead to the expressions for the isothermal sus- 
ceptibility \T and the magnetic specific heat at constant h, Ch, respectively: 

and 

r _ rp ( d 2 G l{NN} \ f da i{NN} \ 

Moreover, for the independent variables h and T, the mixed derivatives of 
the Gibbs potential must satisfy the identity: 



d 2 G if N N\ d 2 G 

{62) 



T i{NN} _ <J (j i{NN} 

dTdh ~ dhdT 
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Taking into account that eqs. (1251) and ([2"g]) are fulfilled, we can write eq. 
(J22J) in the form of Maxwell relation: 



In practice, the validity of relation ([33]) can be checked for a given z, h and T 
numerically only, because the analytical expresions for m and &i{NN} (eqs. ([8]) 
and ( |T9l) . respectively) are too complex. 

The equations (I26l 1281 and [301 EI]) allow us to study the basic properties 
of the model from the free-energy expression, whereas the equilibrium con- 
dition (122]) . and the existing EFT equations in H-K approximation (0 [TT]) . 
remain fulfilled. 



5 Discussion of the results and conclusion 

In the paper a method of derivation of the Gibbs free-energy within H-K ap- 
proximation has been shown. In calculations one has to solve the differential 
equation ([23]) for the variational parameter A. In order to perform a numer- 
ical test for the theory we have calculated the A-parameter vs. temperature 
for h = and z — 4. The testing results are gathered in Fig. 1. For compar- 
ison, the magnetization m and NN correlation function c are also presented, 
being the solutions of eqs. ([H]) and ffTT]) . respectively. The numerical solution 
of differential eq. (123]) can be done by iteration procedure. For this purpose, 
it is convenient to start from the Curie point (when A = 0) and to decrease 
the temperature in small steps, for each step calculating successive m and A 
until the ground state. For T = we obtained the value A = 0.6261. Then, 
having A upon T, as well as m, c, and A from eq. (1141) . the Gibbs free-energy 
can be found from eq. (121 p for arbitrary temperature. 

The results of the Gibbs energy calculations per 1 lattice site (i.e. the 
chemical potential) are also presented in Fig. 1 for z — 4. For comparison, 
the chemical potential in the Molecular Field Approximation (MFA) and in 
the Pair Approximation (PA) is also shown. The calculations for MFA and 
PA are based on the theory given in Ref . [33] . It turned out that the curves 





(33) 
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for EFT in H-K approximation and PA are almost the same (within the ac- 
curacy of the line thickness) for the range of temperatures discussed. This 
is a worth-noticing result, because the PA method is usually considered as 
a more accurate approach, among other things, giving Bethe result for the 
Curie temperature. In fact, Curie temperatures for z = 4 resulting from 
EFT (in H-K approx.) and PA are different and they are k^Tc/J = 0.7725 
and k B T c /J = l/(21n2) = 0.7213, respectively (for comparison, k B T c /J = 1 
in MFA) . Explanation of this difference is due to the different magnetization 
curves (i.e. the /i-dependence of the Gibbs energy) in both models. 

As far as the simple cubic lattice is concerned (z = 6), the Curie tem- 
perature resulting from EFT for S z = ±1/2 is equal to k B T c /J = 1.268. 
This value corresponds to the first Matsudaira's approximation [11], and is 
much better then MFA result k#Tc/J = 1.500. On the other hand, Bethe 
analytical result (equivalent to PA [33J ) is k B Tc/ J = l/(21n(z/(z — 2)), which 
for z = 6 gives the value k B T c /J = 1.233. One of the best results for the 
Curie temperature, k-oTc/J = 1.145, has been obtained by Kikuchi [34] by 
means of cluster variational method in cubic-cell approximation. The result 
of Kikuchi is very close to HTSE Curie temperature k B Tc/J = 1.128 |35j, 
which is probably the most accurate value for s.c. lattice. The Kikuchi 
method is worth-mentioning because it also enables calculation of the en- 
tropy, however, by a more complex combinatorial approach. A disadvantage 
of the Kikuchi method is such that it involves many variational parameters, 
according to which the chemical potential should be minimized. In contrast, 
the advantage of the present EFT theory is that only one kind of parameter 
(magnetization) can be chosen to minimize the Gibbs energy. Contrary to 
PA [33J, or Kikuchi metod [34J, in this formulation of EFT the correlation c 
is not a variational parameter, but serves only as a complementary quantity 
being fully determined by independent variables h and T. This simplicity 
of EFT provides much wider perspectives of application at the expense of 
accuracy. It should be stressed that EFT is superior to MFA, where the 
correlations are totally neglected. 

The method presented in the paper enables complete thermodynamic de- 
scription from one root expresion proposed for the cluster density matrix. 
From the reducibility of this matrix, the EFT equations for magnetization 
and NN correlation functions are re-derived. The expression for the Gibbs 
energy in H-K approximation is obtained more correctly than that presented 
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in Refs.[28, 29J. Let us remind that in the existing EFT approaches to the 
Gibbs energy the eqs. (I2B|1 and (l2"2l were not simultaneously fulfilled. In 
the present approach, all thermodynamic quantities can be calculated as the 
successive derivatives of the Gibbs free-energy, together with the extremum 
condition fl22|) . 

We hope that the method can be employed for more complicated situa- 
tions, for instance, for disordered materials (with site and/or bond disorder) 
as well as for thin films and surfaces. As far as higher spins (S > 1/2) are 
concerned, the problem becomes more complex because higher spin moments 
are then involved in all relevant thermodynamic expressions. However, the 
possible implementations of the method should be published elsewhere. 



6 Appendix: 

A set of coefficients for z = 2, 4 and 6 



We define the functions: 



F n = tanh 


'§{nJ + hj 


— tanh 


f (nJ - 


h) 


G n = tanh 


"f (n J + h) 


+ tanh 


"f (nJ - 


h) 



(Al) 



and 



X n = In {2 cosh [f (n J + h)] } - In {2 cosh [f (nJ - h)] } 
Y n = In {2 cosh [f (nJ + h)]\ + In {2 cosh [f (nJ - h)] } 



The temperature and field dependent coefficients from eqs. (El), (ELL 
are the following: 

z = 2 case: 



(A2) 
and (USD 



a = 


g (^1 + *b) 


ai = 




a 2 = 


^ (^1 - ^0) 



(A3) 
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" 16 
h = -G 1 



(A4) 



c = -(Y 1 + Y ) 

ci = X 1 (A5) 
c 2 = Y 1 -Y 



a = — (F 2 + 4Fi + 3F ) 

a 1 = ^{G 2 + 2G 1 ) 

a 2 = | (F 2 - F ) 
03 = (G 2 - 2Gi) 
a A = \ (F 2 - 4Fi + 3F ) 



(A6) 



60 = 


^(G 2 + 2 Gl 


61= 


I (F 2 + Fx) 


6 2 = 




63 = 




64 = 


^ (G2 " 2d) 



(AT) 
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z = 6 case: 



c = -(Y 2 + 4Y 1 + 3Y ) 
c 1 = i(X 2 + 2X 1 ) 

c 2 = \ (Y 2 - Y ) (A8) 

c 3 = 2(X 2 -2X 1 ) 
c 4 = Y 2 - 4Y 1 + 3F 

a ° = T^o ( F 3 + 6F 2 + 15Fx + lOFo) 
ai = ^(G 3 + 4G 2 + 5G' 1 ) 

a 3 = ^(G 3 -3Gi) (A9) 
a 4 = v ( F 3 - 2F 2 - F 1 + 2F ) 

o 

a 5 = ^(G , 3 -4G' 2 + 5G' 1 ) 

«6 = ^ (F 3 - 6F 2 + 15Fi - lOFo) 



60 = ^(G 3 + 4G' 2 + 5G' 1 ) 

61 = ^ (3F 3 + 8F 2 + 5Fi) 

6 2 = |t(3G' 3 + 4G' 2 -G' 1 ) 
b 3 = ^(F 3 -F 1 ) 

h = ^(SG 3 -AG 2 -G 1 ) 

14 



(A10) 



b 5 = - (3F 3 - 8F 2 + 5Fx) 
6 6 = i(G 3 -4G' 2 + 5G' 1 ) 



c o = A(*3 + 6^ + 15^ + 10^) 

ci = 4 (X3 + 4X2 + 5X0 
lb 

c 2 = ^(Y 3 + 2Y 2 -Y 1 -2Y ) 

03 = ^X3-3x0 (All) 

c4 = ^(y 3 -2y 2 -^i + 2yo) 

c 5 = 3(X 3 -4X 2 + 5X 1 ) 

C 6 = Y 3 -QY 2 + 15Y 1 -10Y 
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Figure 1: Illustration of the variational parameter A/ J, magnetization m, NN 
correlation function c and the chemical potential G/NJ upon the reduced 
temperature h^T j J for z = 4 and h — 0. 
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